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ABSTRACT 

Spectropolarimetry from the near IR to the far UV of light scattered by dust provides a 
valuable diagnostic of the dust composition, grain size distribution and spatial distribution. 
To facilitate the use of this diagnostic, we present detailed calculations of the intensity and 
polarization spectral signature of light scattered by optically thin and optically thick dust in 
various geometries. The polarized light radiative transfer calculations are carried out using 
the adding-doubling method for a plane-parallel slab, and are extended to an optically thick 
sphere by integrating over its surface. The calculations are for the Mathis, Rumple & Nordsieck 
Galactic dust model, and cover the range from 1 to 500 A. 

We find that the wavelength dependence of the scattered light intensity provides a sensitive 
probe of the optical depth of the scattering medium, while the polarization wavelength 
dependence provides a probe of the grain scattering properties, which is practically independent 
of optical depth. We provide a detailed set of predictions, including polarization maps, which 
can be used to probe the properties of dust through imaging spectropolarimetry in the near IR 
to far UV of various Galactic and extragalactic objects. In a following paper we use the codes 
developed here to provide predictions for the dependence of the intensity and polarization on 
grain size distribution and composition. 



Subject headings: dust, extinction — radiative transfer — scattering — polarization — ISM: 
globules — methods: numerical 



1. Introduction 

Light scattering induces polarization, and the intensity and polarization of the scattered light 
tells us the properties of the scattering medium. Continuum radiation can be scattered over large 
volumes of space either by free electrons, or by dust grains in either diffuse or compact clouds. Electron 
scattering is wavelength independent (below X-ray energies) and follows a simple scattering phase function 
(oc 1 -|- cos^ 6scat), but dust scattering is much more complicated since the scattering efficiency and the 
phase function can be strongly dependent on wavelength, grain size, and grain composition. Thus, light 
scattered by dust carries valuable information on the dust properties, and the purpose of this paper is to 
facilitate the extraction of this information. 

Dust reflected light is observed in various Galactic objects, either from large and diffuse reflection 



nebulae (e.g. Witt et al. 1992, 1993; Gordon et al. 1994; Calzetti et al. 1995), or from more compact, but 



still spatially resolved objects, such as dense molecular clumps (e.g. the Eagle nebula, Hester et al. 1996), 
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protoplanetary disks (the Orion nebula, O'Dell, Wen, & Hu 1993; Bally et al. 1998| ), cometary knots in 
planetary nebulae (Balick et al. 199S; Burkert & O'Dell 1998 ), and bipolar nebulae (the Egg nebula, Sahai 
et al. 1998a, b). Dust reflection is also inferred in some compact unresolved objects (e.g. Schulte-Ladbeck 
et al. 1992; Rees & Sitko 1996), and is also likely to be present in young stellar objects with circumstellar 
disks (e.g. Chiang & Goldreich 1997). 

Dust reflected light is also observed in various extragalactic objects. In particular in active galaxies, 
where dust obscuration plays a key role in unification schemes (Antonucci 1993; Urry & Padovani 1995). 
Specific evidence for scattering is present in Seyfert 2 galaxies (e.g. Antonucci & Miller 1985 ; Miller, 



Goodrich & Mathews 1995; Tran 1995a, b; Capetti et al. 1996), radio galaxies (e.g. Jannuzi et al. 1995| ; 



Dey et al. 1996; Cimatti et al. 1998; Tran et al. 1998; Hurt et al. 1999), broad absorption line quasars (e.g. 

Hines et al. 1999), red quasars (e.g. Brotherton 



199E; Goodrich & Miller 1995; Schmidt 



Cohen et al. 

et al. |1998; De Breuck et al. 1995), and ultraluminous IR active galaxies (e.g. Wills et al. 1992; Hines & 



Wills 1993 ; Hines et al. 199£ ). In some of these objects the polarization rises towards shorter wavelengths, 
strongly suggesting optically thin dust scattering. 

Despite its potential usefulness, almost no data currently exists on spectropolarimetry of spatially 
resolved scattered light in Galactic objects (see a rare exception in Code et al. 1996). Fortunately, significant 
amounts of spectropolarimetry, though mostly spatially unresolved, are available for extragalactic objects, 
in particular at the rest frame UV (for high z objects), where most of the spectral polarization features are 
predicted, and which therefore provides the strongest diagnostic power. 

In this paper we provide a detailed set of predictions for the wavelength dependence of the intensity, 
percent polarization, polarization plane, and circular polarization, of light scattered by dust in the near 
IR to the far UV. These properties are explored here for a given dust model over a range of scattering 
geometries for either optically thin or optically thick dust. In a following paper we explore these properties 
for various dust models, as characterized by the grain size distribution and grain composition. These 
predictions, together with new high quality spectropolarimetry in the near IR to the far UV of reflection 
dominated objects, should allow significant constraints on the properties of dust in various Galactic and 
extragalactic environments. 

Earlier calculations of polarized radiative transfer in a dusty optically thick medium were made by 
White (1979a,b), Kartje ( |1995| ), Code & Whitney ( |1995| ), and Wolf & Henning ( |1999| ; see also optically 
thin calculations for radio galaxies in Manzini & di Serego Alighieri 1996). The calculations presented 
here improve on these earlier calculations either by the use of an accurate radiative transfer, rather than a 
Monte Carlo calculation, or by using updated dielectric functions of silicate and graphite grains, by using 
the actual, rather than approximated scattering phase matrix, by extending our results down to 500 A 
(relevant for high redshift extra galactic objects), and by the detailed set of results presented. 

The paper is organized as follows. In §2 we provide a detailed description of the theoretical approach. 
The results for scattering by optically thin dust, and by an optically thick slab and sphere are presented in 
§3, and the paper is summarized in §4. A detailed description of the numerical implementation is provided 
in the Appendix, where we also provide tests of the codes based on a comparison with accurate analytic 
results for electron scattering, and with earlier numerical results for dust scattering. We also show in the 
Appendix a comparison with the much simpler case of an electron scattering sphere. 
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2. Theoretical approach 



The main purpose of this paper is to study in detail the wavelength dependence of scattering and 
polarization of light incident on optically thick dust in simple geometrical configurations. The distance 
from the illumination source to the dusty nebula is assumed here to be much larger than the characteristic 
length corresponding to an optical depth r 1 within the dusty nebula, so one may use the parallel 
illumination assumption. Dusty slabs or spheres are the simplest possible dust configurations to consider, 
and are relevant for the various optically thick scattering objects mentioned in §f . 

Radiative transfer in a plane-parallel slab is a standard problem in the theory of radiative transfer 



(Chandrasekhar |1960| ) . To solve this problem we use the adding-doubling method (Hovenier 1971 ; Hansen 



& Travis 1974; Evans & Stephens 1991), which is efficient and numerically stable up to large optical depths. 



We follow the formulation of the adding-doubling method proposed by Evans & Stephens (1991). To solve 
the case of an optically thick dusty sphere, we assume that the geometrical depth at which t ^ 1 is much 
smaller than the radius of the sphere, so the surface of the sphere can be treated as a collection of a large 
number of small plane-parallel slabs with a horizontal size much less than the sphere radius, but still 
much larger than the vertical depth at which r ~ 1. The scattering properties of the sphere are obtained 
by integrating the local solution for the slabs over the sphere surface. Note that the surface integration 
approach can be applied to practically any optically thick geometry. 

Light scattering by a sphere can also be calculated using the Monte Carlo method (Witt |l977 ), as 
was done by Code & Whitney (|l995|). However, the efficiency of the Monte Carlo method drops steeply at 



large optical depths, typically for r > 5-10 (Fischer, Henning, & Yorke |1994| ; Code & Whitney |1995| ; Wolf 



& Henning 1999), and it does not allow the accuracy provided by an exact numerical radiative transfer 
for a large optical depth. However, the disadvantage of the numerical scheme employed here for a sphere 
is that it applies only for t 1, unlike the Monte Carlo method which allows one to treat spheres with 
intermediate optical depth. 

In the following subsections we provide a detailed description of our theoretical approach. We first 
briefly review the basics of single scattering by ensembles of spherical grains (section ^.l| ), introducing the 
phase matrix, albedo, and other parameters relevant for multiple scattering. In section |2.2| we describe the 
polarized radiative transfer in a slab, which provides the basis for this and forthcoming studies. In section 
we formulate the problem and introduce the main parameters. We next derive the expression for 



2.2.1 



the scattering matrix (section 2.2.2) in a form suitable for efficient implementation of the adding-doubling 



method, which is described in section 2.2.3. In section 2.3, we describe how to calculate the scattering 



characteristics of an optically thick sphere using the solution for a slab. 



2.1. Single scattering by dust grains 

In order to solve the radiative transfer in an extended dusty medium, it is necessary to obtain first 
the scattering properties per unit volume. If the mean distance between dust particles in the medium 
considerably exceeds the typical grain size, and the particles are randomly distributed (as applies for cosmic 
dust), then there are no coherence effects and the particles scatter the light independently. The scattering 
also does not change the frequency. The transfer in larger, optically thick volumes, is studied by solving 
the equations of radiative transfer, using the results for the single scattering problem. In this section we 
introduce the main parameters which describe single scattering; the intensity vector, the phase matrix and 
the single-scattering albedo. We next present specific expressions of these parameters for the cases of Mie 
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and Rayleigh scattering relevant for the present study. 



The intensity and polarization of a beam of light are fully described by the 4-vector I of Stokes 



parameters {I,Q,U,V) (Chandrasekhar 196C; Bohren & Huffman 1983), defined with respect to a 
certain plane of reference. The parameter / describes the total light intensity, Q and U describe the 
linear polarization, frequently expressed through the degree of linear polarization pi ~ \/ + W^/I 
and the position angle of the direction of linear polarization with respect to the local meridian plane 
0p = i arctan([//(5), such that cos20p has the same sign as Q. The last parameter, V, describes the 
circular polarization, and the degree of circular polarization is Pc — V/I. Note that the intensity vector I is 
actually a pseudo-vector, since it does not transform like a true vector under a coordinate transformation. 
It is possible to build the so-called radiation density matrix p in terms of the Stokes parameters (Dolginov, 



Gnedin, & Silant'ev 1995) 



U - 



Q 

iV 



U- 
I 



iV 

Q 



(1) 



which is a tensor in coordinate space. However, we base our calculations on the the intensity vector I, as is 
common for cosmic dust studies. 

The scattering from a small volume is obtained by summing over the scatterings of the grains within 
that volume. We assume a model of spherical grains, which is frequently used in various applications. The 
approach we develop is applicable not only for classical homogeneous spherical grains (Mie particles), but 
also for spherical grains of multilayer or composite structure. 

If a grain of radius a is illuminated by parallel radiation of intensity Ijn at a wavelength A, then the 
intensity of the radiation scattered by the grain in the direction 0scat and measured at a distance r in the 
far-field (r ;» A) is 

(2) 



Iout(Bscat, a, A) = """^ ' 7'(9scat, a, A) li, 



47^7-2 

where crsca(a, A) is the scattering cross-section of the grain and ^^(Oscat, a, A) is a 4-by-4 matrix called the 
phase matrix. 

Since energy can be removed from the incident radiation by scattering as well as by absorption, it 
is important to introduce also the absorption and extinction cross-sections; (Tabs and CToxt = fsca + Cabs, 
respectively. Then we may define the so-called single-scattering albedo uj — (Jsca/ccxt, which gives the 
probability that a photon hitting a grain will be scattered. It is common to use the dimensionless efficiency 
factors for scattering, absorption and extinction, defined as ratios of respective cross-sections to the grain 
geometrical cross-section, equal in our specific case to na^: 



^ <-^sca 

^SCa n" ; 



Qahs — 



Cabs 



^ _ CToxt 
^; cxt — o" J 

na 



(3) 



and thus uj = Qsca/ Qc 



The phase matrix V contains the full information about the angular distribution and state of 
polarization of the scattered radiation. For spherical grains, V has an especially simple expression (Van de 



Hulst |1957| ; Bohren & Huffman |1983D 



P(escat,a,A) 



-Pl(6scat, a, A) /'2(0scat, a. A) 

^'2(6scat, a. A) Fi(9scat, a. A) 











PsiQscat^a, X) P4(6scat, a. A) 

-P4(Bscat,a, A) P3(escat,a, A) 



(4) 
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where only three of the four matrix elements are independent since: = P| + P| + P|. The scattering 
plane, defined by the incoming and outgoing beams of light, is the plane of reference for the vectors of 
the Stokes parameters. The element Pi is the phase function, which describes the probability that the 
unpolarized incident light will be scattered in any specific direction. The phase matrix V is therefore 
normalized such that 

dn 

The angular dependence of the phase function Pi can be characterized by the so-called asymmetry 
parameter g: 

g^icose) - / cos Pi ^, (6) 

where g ranges from -1 (complete backscattering) to 1 (complete forward scattering) with g = 
corresponding to isotropic scattering. 



2.1.1. Mie scattering 



Light scattering by a spherical homogeneous particle of an arbitrary radius a is given by the Mie series 
solution (the formulae presented here are valid also for the case of a concentric multilayer sphere). The 
elements of the phase matrix for a Mie-scattering particle are (Bohren & Huffman 1983 ): 



P - -^(1 ^2 P + I Si |2), P2 ^ ^(1 52 P - I Si n 



^3 = 3^(^2^1+^2^!*), 



2 — 
P4 = 



2i 



■iS^Si-S2Sl), 



(7) 



where x = l-najX is the size parameter, and S'l and ^2 are the scattering amplitudes, in general complex 
functions, for which the Mie theory gives (Bohren & Huffman 1983): 

" 2fc + 1 „ , , A 2fc + 1 



1 ^ ^ k=l 



kik + 1) 



[akTk{y) + hkTikiy)]- 



(8) 



where ^ = cosGscat, the functions nk{p) and Tk{fJ.) contain the information about angular dependence, and 
ak and bk are the so-called scattering coefhcients. The angular functions iTkifJ,) and Tfc(/i) can be calculated 
from the simple recursion relations (Bohren & Huffman 1983): 



2k 



k 



7rfe+2 



-^J■^Tk+l 



-TTfe, Tk+2 ^ {k + 2)/i7rfe+2 - (fc + 3)7rfc+i 



(9) 



k + 1' ^ fc + 1 

starting from 7ro=0 and 7ri=l. The coefficients ak and 6fc, which are generally complex, characterize the 
contribution of various multipole modes to the scattering process. These are dependent on a, A and the 
grain material complex refractive index n ^ rir + irii. The specific expressions for ak and bk, which are also 
calculated with recursion relations, may be found, e.g. in Bohren & Huffman ( 1983| ). 

Once ak and bk are known, the efficiency factors Qsca and Qcxt, a-nd the asymmetry parameter g, are 
given by (Bohren & Huffman 1983): 

2 " 



Qsca = ^^(2fc-f 1)(| ak 

^ k=l 



bk P) 



Qcxt = — 5I(2fc + 1) ^{ak + bk}, 



■E 

k=l 



k{k 



'^Wka*k+i + bkb*k^i} + 



k=l 

2k 



1 



k{k + l) 



"iiiakbl} 



(10) 
(11) 
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Computational implementation of formulae ([zj-pl]) can be found in Bohren & Huffman (1983) and Wiscombe 



(1979) 



2.1.2. Rayleigh scattering 

The Rayleigh scattering approximation applies when the grains are very small, i.e. a; ^ 1, and in 
addition x\n\ <^ 1. The above conditions mean that the particles interact with the radiation mainly as 
electric dipoles with a negligibly small contribution of higher multipoles. The elements of the phase matrix 
in this case can be obtained from those for Mie scattering [eqs (0)] by retaining in the expansions (||) the 
terms responsible for the electric dipole interaction, which gives the very simplified expressions (Bohren & 
Huffman |198|): 

Pi = 1(1 + C0S2 e^cat), P2 = -| Sin^ Gscat, ,-^2) 

The phase matrix V is a function of the scattering angle 6scat only, unlike the phase matrix for Mie 
scattering particles which depends also on grain properties: a, n, and A. The expressions for the optical 
efficiencies also become simpler: 



3 



1 



Qext = 4x Im <; } ■- (13) 



where Im stands for the imaginary part. Note that for Rayleigh scattering by free electrons (Thomson 
scattering) the elements of the phase matrix are the same as in (|l2|). 



2.1.3. Single scattering by a mixture of size- distributed grains 

In the previous sections we considered light scattering by a single grain. A dust mixture contains 
grains made up of various species and a range of sizes. Let /'(a) be the size distribution function for the 
ith dust species, that is the number density of grains of type i having sizes within the interval a to a+da. 
The scattering parameters of a volume element is obtained by averaging the single grain scattering over the 
dust species and the grain-size distributions, which yields the following expressions: 

p.p. _ EJo°°^«' QLa(a.A) P,'(9scat,a,A) Pja) da 

Jo QscAa, A) P{a) da 

n _ E»/o°°7ra2 QlcM,>^) fja) da 

Li Jo 7ra^r(a)da 

n M^ _ E^/o°°^"^ Qcxt(a,A) f{a) da 
Li Jo TTfl^ r(a) da 

. _ E» ir ™^ Qsca(a, A) g(a. A) /'(a) da 

L.r™^QU«,A)f(a)da ' ^'^ 



EJo°°^^' QLaKA) r{a) da 
LJo°°™' Qoxt(«,A) fia) da' 



(18) 
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where the sum index i runs over all grain species. It is useful to define the scattering and extinction 
coefficients: 

/>oo 

Asca(A) = V / TTa' QL,{a, A) f{a) da, (19) 

Acxt(A) = V / ™' gLt(«, A) /'(a) da, (20) 

measured in units of length"^. The inverse quantities k~J^ and fc^J^^ represent the mean free path of photons 
for scattering and extinction. 

2.2. Scattering by an optically thick dusty slab 

In this section we formulate the solution for polarized radiative transfer in a a plane-parallel dusty slab. 
Since dust scattering does not change A, we omit the explicit dependence of the intensity and polarization 
on A. 

2.2.1. Mathematical formulation of the problem 

Let a dusty slab of optical thickness tq be illuminated by unpolarized radiation incident from a direction 
6*0 with a flux (see Fig. |l|). The equation of polarized radiative transfer for the slab is (Chandrasckhar 

ileol ): 

^ = -I(r, /i, *) + F(t, /i, I) + S(t, fi, $), (21) 
F(t, /i, $, I) ^ J Tifi, <&, m', <f ') I(t, /i', $') V da>', (22) 

S(t, M, *) = £ ^0 $, Mo, I'o) 1 + {I - ^) B{T) 1, (23) 

where I is the intensity vector of the diffuse radiation at an optical depth r, going along the direction (0, 
$) or (/i, $); where /i is the cosine of zenith angle 0, measured from the upper normal for the upward 
intensity and from the downward normal for the downward intensity (see Fig. |l|); $ is the azimuthal angle, 
measured counterclockwise from the projection of the incoming central source beam onto the slab. Note 
that all scattering quantities depend on the azimuth difference rather than on the azimuth itself. For the 
incoming beam $o = and /iq — cos 6q . 

F is the scattering integral, expressing the angular transformation of the diffuse radiation due to 
scattering. T is the 4-by-4 scattering matrix that relates the intensities of incoming and outgoing radiation, 
taken in their specific reference planes (see section 2.2.2| for more details) . 



S is the source function 4-vector, which includes a 'pseudo-source' of single scattered radiation of the 
central source and unpolarized thermal emission of dust. The latter may be approximated by a black-body 
radiation B of temperature T. 1 is the unity 4-vector (1, 0, 0, 0). Note that in this study we are interested 
in light scattering in the UV and the optical, where the contribution of the dust thermal emission is 
negligibly small. We therefore exclude this effect below. 
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To complete the formulation of the problem, we should specify the boundary conditions for the intensity 
on the top and bottom faces of the slab (Fig. Q): 



1(0, 



(24) 



where Ij (ij,) is the intensity of diffuse radiation from external sources incident on the top (bottom) of the 
slab. For a slab with a large optical thickness tq ^ 1, instead of the second boundary condition in (^J) we 
expect that the intensity I(r, ^, $) at r — > oo be limited. Note that in this study the intensities and ij^ of 
external diffuse radiation in the optical/UV are expected to be negligibly small compared to the intensity of 
the central source, so we may set them to zero. However, we retain I^f and ij, until the end of this section 
for purpose of generality. 

In the following sections we describe how the problem of radiative transfer in a dusty slab can be solved 



efhcicntly by the doubling and adding method (Hovenier 1971; Hansen & Travis 1974) 



2.2.2. Scattering matrix 

The solution of the radiative transfer problem with the doubling and adding method requires a suitable 
presentation of the scattering matrix in terms of angular dependence. Fig. |^ displays the geometry of 
scattering for any particular point inside the slab. We have an incoming (^in) and outgoing (/J.out, 'i'out) 
beams with the scattering angle between them Oscat- The theory of single light scattering provides the phase 
matrix .7^ as a function of scattering angle [e.g., see (||), (0), ( [l^ ) and (^], where here the scattering plane 
is the natural plane of reference. However, both incoming and outgoing beams have their own reference 
planes, defined by the Z-axis and a respective beam. Thus, we should perform a few transformations in 
order to get the scattering matrix J- from the phase matrix V for any particular beam (Chandrasekhar 
1960| ): 

•^(Mout, *out, Aiin) T^{i2 - Tt) V{cOS Oscat) 7^(«i), (25) 

where TZ is the polarization rotation matrix, which is a function of the rotation angle i: 



n(i) 



1 








" 





cos 2i 


— sin 2i 








sin 2i 


cos 2i 














1 



(26) 



Note that this form of TZ is valid for the intensity vector representation (/, Q, U, V), adopted by us. In ([25|), 
the matrix TZ first transforms the intensity vector of the incident beam from its reference plane A™OZ 
into the scattering plane A™OA°^^ by rotating through the angle ii between the two planes (see Fig. H]). 
The matrix V then gives the intensity of the emergent beam, scattered by Gscat- Finally, the matrix TZ 
transforms the intensity of the emergent beam from the scattering plane into the reference plane of the 
beam ZOA°^^ by rotating through the angle 12 — tt. Note that the transformations (|2^) affect the linear 
polarization parameters Q and U only, leaving the mean intensity / and circular polarization parameter 
V unchanged. Only the position angle 9p of the polarization plane is changed, while the degree of linear 
polarization p\ = \J + jl is unchanged. The scattering angle Oscat, and the angles i\ and ii are 
calculated from the angles /iin, /iout and ^out: 



cos ©s 



MinMout + \/(l -Mhi)(l -Mout) 



COS $n 



(27) 
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COSZi 



C0SZ2 



A^out x/l-MS"-Min\/l - Mout COS ^out 
sin Oscat 

A^in \/l - Mout - A^out y/l - fJ-l, COS ^out 
sin Oscat 



(28) 
(29) 



For computational purposes, it is useful to expand the scattering matrix in a Fourier series in azimuth 
angle $out: 



^(Atout, $out, Min) = ^ [■?Vi(A*out, Mill) COS m$out + (Atout , Min) Siu m$o 



(30) 



Then the Fourier modes J-^ and J-^ can be calculated with a computationally efficient fast Fourier 
transform (FFT) for each pair of jj-m and A^out: 



•^m(A*out J A'in) + ^•^mlA'outj A^in) 



(2 — S„m) 

N 



N-l 
n=0 



^(A*out, '&outJ Min), 



(31) 



where = 2TTn/N for n ranging from to where N must be a power of 2. 



In this paper we deal with mixtures of spherical grains for which the phase matrix V has a simple form 
(^. For this type of V, the scattering matrix T has the following symmetries: 1. Changing the sign of 
/Ltout and /Ltin results in a change of sign of the off-diagonal 2-by-2 blocks of the matrix. Thus, only half the 
number of /is is needed to compute the whole scattering matrix. 2. The upper left and lower right 2-by-2 
blocks of are even functions of <&out, while the upper right and lower left blocks are odd functions. Thus, 
the scattering matrix for tt < $out < 27r can be computed from the values for < $out < tt. 

From the last symmetry it follows that the cosine matrices (sine matrices J-^) have off-diagonal 
(diagonal) blocks of zeros. This means that we can compose a single 4-by-4 scattering matrix for each 
azimuth mode and pair of zenith angles. Defining a Fourier intensity vector of the form I = {I'^, Q'^ , U^, V^), 
containing the cosine azimuth modes of the / and Q Stokes parameters and the sine modes of U and V, 
then the Fourier scattering matrix can be represented through the Fourier terms from the FFT as 



^ TM (A^out ; A^in) 



•r-c 
11 


-r-c 
12 


-r-s 
•^13 


-r-s -1 
■r 14 


•r-c 
21 


-r-c 
22 


-r-s 
•^23 


-r-s 
24 


-r-s 
-^31 


•r-s 
■^32 


-r-c 
•^33 


-r-c 

34 


-r-s 
J- 41 


_ •r-s 
•^42 


-r-c 
•^43 


-r-c 

^44 J 



(32) 



This form of scattering matrix is especially suitable for efficient implementation of the doubling and adding 
method (see next subsection). 



2.2.3. The doubling and adding method 

Fortunately, there is an approach to solve the radiative transfer problem (|2^-|2^), which does not 
require a direct solution of the complicated integro-differential equation (|l]). This approach is based on 
the general principle of the linearity of interaction of electromagnetic radiation with a scattering medium. 
Thus, the intensity of radiation leaving a slab is linearly related to the intensity of the radiation incident on 
the slab: 

ji ^-^iiji_^-^iTiT^S^ ll=T^^ll+n^^\i+S\ (33) 
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where Ij (If, ) are the intensity vectors on the top (bottom) of the slab, T^^ and T^^ (TZ^^ and 7^^^) are 
4-by-4 transmission (reflection) matrices, and S''-'^ are the source vectors. The arrows denote the downward 
(I) and upward (|) directions of the beams. In ( |3^ ) we use a short notation, where the product J=AI of a 
matrix A and a vector I means: 



J»(a*,*) = 



1 



fi' d^' d$'. 



(34) 



where the index i, indicating the Stokes parameter, runs form 1 to 4. A similar rule is valid also for a 
matrix by matrix multiplication, e.g. for C^AB 



Cii'(^,/iin, $) = 



1 



n 1 n2TT 

'o Jo 



A^^" (a*, m', * - *') iP-', Aiin, $') 



U" = l 



^i' d/i' d$'. 



(35) 



Equations ( |3^ ) can be transformed into a form more suitable for computations by the following steps: 1) 
expand the matrices and vectors in (|3^) in Fourier series in azimuth $ and choose the optimum Fourier basis 
that would result in decoupling of azimuth modes; 2) discretize zenith angles with a numerical quadrature; 
3) rearrange the respective vectors and matrices to their normal forms as one- and two-dimensional objects. 
Below we describe how these steps are realized for relations (^J) and ( |35| ) . 

Expanding the quantities in (p^-05|) in Fourier series in $ and choosing a representation for the Fourier 
components of vectors as (P ,V^) and for matrices the one like ( ^2| ) (see section 2.2.2), decouples 
the azimuth modes such that for each mode m: 5™=^™!™ and C^"—A'"B™ where 



jr(M) 



0(Ai,Min) 



i' = l 
■ 4 



/i' d/i'. 



(36) 



(37) 



We now discretize the zenith angles using a numerical quadrature, for example, the Gaussian formula. 
The expressions are then transformed into 



9 \ ^ \ ^ Am -rm 

/ . / . "i' ^i' ■^ii' ,jj' i' ,j' ' 
= l i' = l 



(38) 
(39) 



Cii'jj' — 2 ^ ^ fij" Wj" Aii„^jj„ Bi„i, .j„j,, 
j"=i t"=i 

where the indices j and j' indicate the quadrature angles /ij and weights Wj, and run from 1 to n^, the 
number of quadrature angles. 

Finally, it is useful to rearrange the intensity vector I™ for the mth azimuth mode from a two- 
dimensional array into a one-dimensional vector J™, where k = i + 4(j — 1) is running from 1 to 4np. 
If we then respectively redefine the matrices as, say, A™f^, = 2^jiWjiA™, jj, with k — i + 4(j — 1) and 
k' = i' + 4(j' — 1) instead of A^, jj,, the expressions for multiplications of vectors and matrices acquire their 
natural form: 



jrn ^ ^ ^ ^ 



kk" ^k"k'^ 



(40) 



k' = l 



k" = l 
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for J™ = j\rnjm ^^^^^ (jm _ j^rn j^m ^ respectively, where the vectors have a dimension n — AUf, 



We may now return to equations (^) which relates the intensities of the radiations incident onto and 
emergent from the slab. With the new vector basis and new notation, we may rewrite equations ( p3| ) for 
each azimuth mode (index m is omitted for better clarity) as ordinary vector-matrix relations: 



(41) 



Equations ( ^l]) are written for the whole slab. Thus, the matrices T and R and sources S for the whole slab 
provide the solution of our radiative transfer problem. For an optically thin layer Tin ^ 1, in which single 
scattering prevails, we may write explicit expressions for the reflection and transmission matrices and the 



source vector (e.g. Evans & Stephens 1991) 



WTin 1 + do„ 



4 

1 + Son 



I \Tk'^ Su'S, 



Tin 
Mi 



Ifcfc'— Sii'Sjji — 



Sii' Sjj' 



5ii'6jj' 



~ 1 + <5or, 



4 



I J'mi Mi? Mi') 



(42) 

(43) 

(44) 

(45) 
(46) 



Mi Mi 
where m is the azimuth mode, are the Stokes parameter indices, are the quadrature angle indices, 
k ~ i + 4{j — 1) and k' = i' + 4(j' — 1) and | S |™ is the mth Fourier mode of the source term S [eq. (|2^)] of 
the radiative transfer equation. 

The next step is to calculate the quantities f, R and 5 for the whole slab starting from those of the 
optically thin layer. Applying the relations (^) for two adjacent layers and eliminating the intensities at 
the common boundary results in the expressions for T, R and S for the layer combined from the two layers. 
This is the essence of the adding method. It thus allows us to get the properties of the combined layer 
through the properties of the top (T) and the bottom (±) layers: 



Ri^ = R^^ +T\}M^^Ri^TV 



(47) 



where M^^'^^ are the matrices which represent multiple reflection factors. Physically, the formulae ( ^7|) 



may be interpreted through the concept of multiple reflected beams (Hansen & Travis 1974 ). The special 
case of the adding method, when formulae ( p7[ ) are applied to two identical layers with the same optical 
thickness t and dust composition (albedo a), and scattering matrix J-"), is called the doubling method. 

Thus, in order to calculate the matrices R, T and source vectors 5* for a general, inhomogeneous, and 
optically thick slab, one needs to: 1. Divide the slab to thinner layers, which could be approximated as 
homogeneous; the number and thicknesses of the layers are chosen, depending on the required accuracy. 2. 
For each homogeneous layer, apply the doubling method, that is, first calculate the quantities R, T and 
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S for a sufficiently optically thin sublayer, normally ri„ = 10^^ — lO""'', with the formulae ( |4^-|46| ). Then 
build up the quantities R, T and S for the whole layer in a doubling manner using equations ([47|). 3. Build 
up the quantities R, T and S for the whole slab by reapplying the adding formulae (^ ) for all the layers of 
the slab. 



Note that the pseudo-source of the single-scattered radiation from the central source [equation (p3[)] 
is exponentially dependent on the optical depth, thus violating the homogeneity of the layer, required for 
doubling to work. In this case, however, it is possible to use Wiscombe's ( 1976 ) extension of the doubling 
method, incorporating this kind of sources. 



Once the matrices R, T and sources S are computed, we may calculate the Fourier modes of the 
intensities of the radiation emergent from the slab, j-^ and ll^j^ '^^^h formulae (^l|). To get the 
intensities themselves, we should convert the Fourier modes back into the azimuth space: 



rd>^-V"'" T^i" J cosm$, i^l,2\ 
rT _ Y-"'" ft™ I' cosm$, i=1.2 



(48) 



where i is the Stokes parameter index, j is the quadrature zenith angle index and rim is the number of 
azimuth modes. 

Using relations (|4^), we may derive expressions for the intensity vectors /-'' and 7^ at any depth inside 
the slab 

ji ^ j^fll ^1 _^ fU jl j^iT ^_ fjT/T^ , /T ^ mTT + rj^jT _^ jiU ^gl^ ^ T|^//) , (49) 

where the matrices T and R and sources S correspond to the layers above (T) and below (_L) the level, and 
the matrices M are defined through matrices R in (^7|). The Fourier modes of the intensities and D can 
be converted into the azimuth space by means of the formulae like (^) . 

The polarization properties of the radiation reflected by the slab are represented by the degree of linear 
polarization, p\^^, the position angle, 9p, and the degree of circular polarization, pjij^, 

pL-^^^^^, ^*=^arctan^, p^, = |, (50) 

where lj= {It,Qt,Ut,Vt), cos26'p should have the same sign as Qt, and all the quantities in ( |50| ) are written 
for the radiation, reflected from the top. Similar expressions may be written for the transmitted radiation, 
escaping from the bottom. 

Other useful characteristics of the radiation leaving the slab are the fluxes of the radiation emergent 
from top and bottom: 

= /o 7o ItV] (A*' *) = 2n jy;^, w,^^, //°^., , z = 1, . . . , 4, 

where only zero-th azimuth mode of respective intensity need be used. Only the Stokes parameters / and 
Q (i=l,2) of the fluxes fJj^j and F^^^ can be non-zero. The fraction of the incident flux which is reflected is 
called the plane albedo: 

AM - (52) 
where -fo is the flux of the incident unpolarized radiation of the central source, coming at zenith angle fio- 
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2.3. Scattering by an optically thick dusty sphere 

In this section we apply the dusty slab calculations to calculate the scattering from a dusty sphere. The 
sphere is assumed to have r ;2> 1, so that the thickness of the t ^ 1 surface layer is much smaller than the 
radius of the sphere, and thus every surface element can be treated with the "local slab" approximation. As 
in the slab calculations, we assume that the illumination source distance is much larger than the radius of 
the sphere, so the incident rays are assumed to be parallel. The geometry of the scattering is schematically 
represented at Fig. ^. The radius of the sphere is Rs- The flux of the incoming unpolarized radiation from 
the central source is Fq. The coordinate system is chosen in such a way that the X axis is directed toward 
the observer, and the direction of the beam incident from the source is parallel to the XY plane. The 
position of any point S on the sphere surface is then defined by two angles: zenith angle 6 measured from 
the Z axis and azimuth $ measured from the X axis. The angle between the directions from the source and 
to the observer, Oobsj is the same at any point on the sphere surface, and it is equal to the scattering angle 
for single scattering. 

The angular parameters of the incident and emergent beams, are shown at Fig. ^ 

cos6iin(6', $, Gobs) = Min = " smO cos($ + Gobs), (53) 
cos6iout(6',^') = Mout = sin 6* cos$, (54) 

^ //I ^ N cos 9obs + MinAiout ,rr\ 

COS$out(0,$,eob.) = ^=======. (55) 

The intensity vector Igph of the radiation scattered by the sphere is calculated by integrating the local 
intensity 1'°^^ over the illuminated part of the sphere, as seen from the observer direction: 



:,ph= / /^l'°^Moutd5, (56) 



where 1^°^^ is the local intensity of the emergent radiation, as measured in the observer reference plane, 
chosen here to be the scattering plane A™SA°^^. 1^°"^ is derived from the local emergent intensity, ij [see 
equation (^8|)], defined with respect to the local (meridian) reference plane, by rotating it with the matrix 
TZ through a respective angle "0: 

I'°^=7^(v)I^. (57) 



Substituting into (|5^) and ( |57| ) the expressions for the surface element dS = R^ sin 9 d9 d$ and /iout from 
(^), we obtain 

I«ph(i?s, Oobs) = -Rs / ' cos $ d$ / sin^ e d0 ■ I'°^(0, $, Gobs), (58) 

ii°^(0, $, Oobs) = nmo, $)] • ilie, Oobs), (59) 

where the integration limits correspond to the illuminated part of the sphere, as seen by the observer. The 
intensity ij depends on 9, <& and 6obs through ^in, /iout and <i>out (eqs p3|-p5|). V is the rotation angle 
defined through its cosine and sine: 

sin 9 sin $ — cos 9 

cosip =—===, sm-0 = — ^== (60) 

V 1 ^ Mout V 1 ^ M5ut 
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It is simple to show that from the specific azimuthal dependences of various Stokes parameters (see eqs p^ ) 
one gets symmetry relations for the intensities at the surface points, located symmetrically with respect to 
the equatorial plane XOY: 



ItV](^,$,eobs) = l!H(^-e,$,eobs), ^ = l,2, 
lj[^](0,<i>,eobs) = -lt\,](7r-0,$,eob.), * = 3,4. 



Note that these relations are also valid for the intensities 1'°^^. Thus, using relations (^ij) for the rotated 
intensities 1'°'^, we may rewrite the expression ( p8| ) for each Stokes parameter of Igph as follows: 

/sph(i?s,eob.) = 2i?2 cos$d$ /'sin2 0d0-l|°f(0,$,eobs), (62) 

•^f-eob= Jo 

QsphiRs, Qohs) = 2i?2 / ' cos $ d$ / 'sin^ e d0 ■ I|°f (0, $, Oobs), (63) 

C/sph(i?s, Qohs) = 14ph(i?s, Oobs) = 0. (64) 

These symmetry relations imply that the total circular polarization of light scattered by a sphere is zero 
(just as there is no circular polarization from a single scattering by a spherical grain), and the linear 
polarization is 

Psph(6obs) = (65) 

The ratio of the total flux scattered by the sphere to the incident flux {Fq ■ ttR'^) is the so-called spherical 
albedo: 

»i 



^sph = 2 / A{iia) A*o d/io, (66) 

JQ 

where A{iiq) is the plane albedo defined by means of eq. (Ik 



Once the local intensity vector 1'°'^ (^loc, Qioc, f^ioc, Moc) is calculated as a function of the surface 
coordinates (9, $), one can build maps of both linear and circular polarization. For this purpose we define 
the following parameters: the local degree of linear polarization p\°'^, the local position angle 6^°'^, and the 
local degree of circular polarization pjf^ in a usual manner: 



Pl (t^, Wobsj = 7 — fa rf. N ' 



/ioc(0,*,eobs) 

£/ioc(g,$,e 
2 Qioc(0,$,eobs) 



9'-{e, eob.) = I arctan "i'^'f^H^'i , (68) 



r^°<=(f) (T) (=> ^ - Moc(6',$,eobs) , . 

^'''^'''°"^-/ioc(^,$,e.bs)- ^''^ 



3. Results 

This section describes extensive numerical results obtained with the codes PRT and PRT-SPH, which 
are described in detail in the Appendix. The objects we model here are optically thick dusty slabs and 
spheres, composed of a simple mixture of bare spherical graphite and silicate grains having a power-law 
size distribution /(a) ^ a^^-^ for radii 0.005-0.25 /im. This model was introduced by Mathis, Rumpl, 
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& Nordsieck (1977; hereafter MRN) for interstellar medium dust and further developed by Draine & Lee 
(1984). It served as a standard model for interstellar dust for many years, though recent new observational 
data suggest that other dust models may be more viable (see Zubko 1999 for a short review). However, a 
dust model which is consistent with all the new data is still not established, and we therefore chose to make 
our first exploration of the effects of optical depth and geometrical configuration using the standard MRN 
model. 



3.1. Single scattering by MRN dust 

Figure ^ displays the single scattering albedo and the phase function asymmetry parameter, g, of MRN 
dust for the wavelength range 0.05-1 /im, used throughout this section. The optical constants of graphite 
and 'astronomical' silicate were taken from Laor & Draine ( 1993| ). The albedo is typically 0.4-0.6, with 



a significant drop at A < 2500 A, whereas g generally increases towards smaller A. The scattering phase 
function becomes significantly more forward peaked at smaller A due to the increasing size parameter. 
Figure ^ demonstrates that the extinction curve of MRN dust reproduces quite well the mean Galactic 
extinction curve. 

Figure ^ summarizes the scattering characteristics of optically thin MRN dust as a function of A 
and as a function of 0scat- The scattered spectrum is bluer than the incident spectrum at all Ggcat for 
A > 2500 A since the scattering cross section increases with decreasing A. The "bluening" of the scattered 
flux is most prominent at small 6scat since the scattering gets more forward peaked with decreasing A (Fig. 
4). The bluening effect stops below 2500 A for 0scat > 30° due to the sudden drop in the albedo (Fig. 
4), compounded by the rising g. A rather strong scattering peak due to graphite grains is predicted at 
A ~ 2200 A for large Oscat- Another similar, but much broader peak is expected in the far UV (700-1000 A) 
due to both graphite and silicate grains. 

The polarization spectrum p{\) displays a more complicated wavelength dependence than the intensity, 
and it is also strongly dependent on 6scat- Overall, p(A) tend to decreases from 1 /zm to ~ 2000 A, and 
then tends to sharply increase to the far UV. The polarization plane is generally at right angle to the 
scattering plane (defined by the incoming and outgoing rays), excluding a small range of angles at large 
©scat (150° — 170°), where the polarization rotates by 90°, but only over a small range of A, and with a 
very small amplitude. The feature around A=0.14 /im is also predicted for the extinction curve, but it is 
not observed (Fig. ^). Thus, this feature is most likely an artifact due to imperfect optical constants of 
silicate, as discussed by Kim & Martin (1995). 

The angular dependence of the scattered intensity demonstrates how the scattering becomes strongly 
forward peaked with decreasing A. At the longest wavelengths the phase function tends towards the 
Rayleigh scattering phase function limit. The angular dependence of p always shows a clear peak at 
0scat=8O-lOO°. The shape of this peak is not always symmetric, and as mentioned above, p rotates by 90° 
(indicated by negative values in the plot) for a small range of A at large ©scat- 



3.2. Optically thick dusty slab 

The scattering characteristics for an optically thick dusty slab depend on three angles: ^inCA^ii 
9out(Mout) and <i>out (see Figure ||), rather than just Oscat as in the optically thin case. 
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Figure 1^ shows the mtensity It as a function of A for various values of /iin, /iout, and $out- At a given 
^out, Mill, and A, It decreases with increasing /lout, i-e. the surface brightness of the slab decreases towards 
face on view. This effect is largest when is smallest, since then the radiation is incident close to edge 
on, penetrates only a thin layer at the top of the slab, and this layer is effectively optically thin, and of 
low surface brightness when viewed close to face on. The dependence on A is generally weak longward of 
2200 A, with a significant drop at A <2200 A due to the albedo drop. 

Figure || shows P\iJ^)^) for the same parameters. As shown later, the main physical parameter which 
determines the shape of p\iJ^X) is 0scat, which is now a function of /iin, Mout, and <I>out- For example, when 
/iin is very close to edge on (0.0672), and <i>out = 90°, then all values of /iout correspond to 6scat ~ 90°, and 
all the Piin(A) curves in this subpanel nearly overlap. 

Figure ^ shows the position angle 0* for the same parameters. The position angle is almost always 
fixed at right angle to the plane defined by the incoming and outgoing rays. The plane used to define 0* in 
the optically thick case is defined by the slab normal and the outgoing ray, and thus 0* simply represents 
the angle between this plane and the incident/outgoing rays plane. However, this figure is shown since 

jumps by 90° for certain combinations of parameters. This effect is similar to the rotation seen in the 
optically thin case, and occurs here also for large values of Oscat- 

Figure ^ shows the circular polarization pjijj. for the same parameters. The amplitude of p\^^ is 
generally well below 1%, and it can be both negative and positive. When $out=0° or 180°, p\]^—0 for any 
/Xin and /iout- 

The doubling method does not allow one to separate directly the contribution of single and multiple 
scattering to the total polarization. To explore this point we plot in Figure ^pfinl-^) for complete 
backscattering (Oscat = 180°), which occurs for /iin=/iout and $out = 180°. In this case single scattering 
produces no polarization, and the observed polarization is due to multiple scattering only. The typical 
values of Piin(A) are 1-4%, and the polarization plane is generally parallel to the slab normal (see right 
panel of Figure ||). Figure || indicates that Piin(A) is typically 10-20 times larger than obtained for multiple 
scattering, and thus it must be dominated by the single scattering contribution. The low multiple scattering 
polarization is explained by the fact that the radiation field incident on a grain from multiple scattering 
is much more isotropic than the primary source illumination. Most of the photons incident on the grain 
following multiple scattering move roughly parallel to the slab surface, and thus they produce a polarization 
which is generally parallel to the slab normal. 

Figure 02 displays the slab albedo A{\). Its wavelength dependence is very similar to the MRN dust 
albedo (Fig. 0), but the amplitude is significantly lower. At very small /iin the scattering layer is optically 
thin for most /iout: and most photons that scatter upwards escape, while most photons scattered downward 
are absorbed, yielding a slab albedo which is very close to one half the dust albedo. As /iin increases the 
scattering layer becomes optically thicker for all /iout, and the escaping fraction decreases. This effect 
becomes most prominent at the shortest wavelengths where the scattering becomes most forward peaked, 
which decreases the probability the photons will be scattered back and escape. Thus, the typical dusty slab 
albedo falls from ~ 10% in the optical to only a few % in the far UV. 
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3.3. Optically thick dusty sphere 

Figures p"3|-p^ present maps of the linear and circular polarization at A=0.05, 0.2124 and 1 fim for a 
range of Gobs- The illuminated part of the sphere extends to the terminator at <&tcrm = § ^ ©obsi which 
corresponds to a fractional projected illuminated area of (5(0obs) = -^(1 — cosOobs)- 

The polarization at each point on the surface is generally largest at 0obs — 90°, and at most angles 
it is aligned at right angle to the scattering plane (i.e. along the axis of the sphere), indicating that it is 
dominated by single scattering. When Gobs =180° the polarization tends to be aligned radially, indicating it 
originates in multiple scattering which produces a polarization aligned along the local surface normal (see 
previous section). The integrated polarization at 0obs=18O° is obviously zero, as expected from symmetry. 
A similar effect occurs at very small Gobs: where the single scattering polarization amplitude becomes very 
small, and the polarization can be dominated by multiple scattering. For example, when Qohs—'^° and 
A=l fim the polarization pattern is radial, and thus dominated by multiple scattering. But, at A=500 A 
the polarization pattern remains parallel due to the lower dust albedo, and increased forward scattering, 
both of which suppress the effect of multiple scattering. 

The circular polarization maps typically consist of four regions with two boundaries across which 
the polarization changes sign, one is the equatorial line, 0h=O, and the other is the meridional line 



= ^{t^ — ©obs) (see Figures 13-151). There is a critical angle 9^;"^ ~ 60° - 80° at which the circular 
polarization changes sign at all positions, and there is also a critical wavelength A'^"* ~700 A, at 
which a similar sign change occurs. The local amplitude of pc reaches a maximum of around 1-2 % at 
0obs=llO°-13O°. At 6obs=0° or 180°, Pc^O throughout the whole illuminated disk, as expected from 
symmetry. 

Figure |l^ shows the integrated scattering characteristics of an optically thick dusty sphere as a function 
of A and Gobs- This figure is analogous to Fig. ^ which shows the scattering characteristics of optically thin 
dust. The increase in /sph(A, 9obs) with increasing Oobs is due to the increasing illuminated fraction of the 
sphere. The increase in /sph(A, Oobs) is smallest at the shortest A because of the counteracting effect of the 
scattering phase function, which enhances the scattering at small Oobs- Note that the feature related to the 
2200 A extinction bump changes its nature from an emission peak at A ~ 2500 A at large Oobs, to a weak 
absorption trough at A ^ 2100 A at small Oobs- The wavelength dependence of Psph looks quite similar to 
the one for optically thin dust, as further discussed below. 

The spherical albedo as a function of A is displayed in Figure |l^. It is generally rather low, and drops 
from ^ 11% in the near IR to ~ 4% in the far UV. The spherical albedo is very similar to the albedo of a 
slab illuminated at /Zin = 0.5, which is the mean value of /iin for a sphere. 

Figure ^ provides a direct comparison of scattering from optically thin dust and from an optically 
thick sphere, where in both cases the observed spectrum depends on a single angle only (Oscat, or Oobs)- 
The spectrum of the scattered light is quite different in the two cases. It is significantly bluened in the 
optically thin case at A > 2200 A, and is generally grey scattered in the optically thick case. Optically thin 
scattering produces a broad deep centered at A ~ 1500 A, while the optically thick dust scattering shows 
an overall drop at A < 2200 A, with only a mild rise towards 1000 A. In sharp contrast to the intensity 
dependence on optical depth, p{X) is remarkably similar for both cases. The shape of p(A) is almost identical 
for both cases at a given O, and only the amplitudes of p(A) differ by ~ 10%. This similarity reflects the 
overall small contribution of multiple scattering to the polarization. The O dependence of p{X) is also very 
similar for the optically thin and thick cases. Thus, /(A) serves as a sensitive probe of the optical depth 
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of the scattering dust, while p{X) serves as an optical depth independent probe of the dust properties and 
scattering geometry. 

4. Conclusions 

We present extensive numerical results on the scattering and polarization properties of dust in 
various configurations. The solution of the polarized radiative transfer problem is carried out using the 
adding-doubling method, which applies to a plane parallel configuration. This method can be extended 
to other optically thick geometries by an appropriate integration over their surface. The numerical 
implementation is made with a set of codes which take as input a specified set of dielectric functions, 
calculates the scattering phase matrix for a spherical grain of a given size and composition, and integrates 
the matrix elements over a specified dust model to obtain the dust scattering phase matrix. The dust 
scattering phase matrix is the input for the numerical code which implements the doubling- adding method, 
and outputs the Stokes vectors for a large grid of incident and outgoing rays. This grid provides the input 
for the final code which integrates the slab results over the surface of a sphere and outputs the Stokes vector 
as a function of scattering angle from the sphere. 

The accuracy of the radiative transfer calculations was verified by comparison to analytic results 
available for a Rayleigh scattering atmosphere, and by comparison to earlier calculations for an MRN dusty 
sphere. We also used our code to obtain a simple analytic approximation for the flux transmitted through 
an electron scattering slab, which is accurate to 4% at all r. 

In this paper we used this set of codes to explore the effects of optical depth, scattering angles, and 
dust configuration (slab or sphere), on the scattering and polarization properties of MRN dust from the 
near IR (1 iim) to the far UV (500 A). We find that the scattered spectrum provides a sensitive probe for 
the optical depth of the scattering medium, while the polarization spectrum provides a probe of the dust 
model and scattering geometry, which is practically independ(uit of optical depth. We provide maps of the 
linear and circular polarization for scattering from a sphere which may be useful for interpreting imaging 
spectropolarimetry of scattered light in a variety of Galactic and extragalactic objects from the near IR to 
the far UV. 

The MRN model assumes spherical grains, but the fact that extinction in the Galaxy induces 
polarization clearly indicates that this assumption is not valid. How does the grains non sphericity and 
alignment affects the scattering polarization? An accurate answer requires detailed calculations which 
are beyond the scope of this paper. However, we note that the amplitude of polarization induced by 
transmission through an optical depth of unity (< 3%) is about an order of magnitude smaller than the 
polarization induced by scattering. Thus, it is plausible to assume that the observed level of non sphericity 
and alignment will not have a major effect on the results presented here. 

In a following paper we explore the polarization signature of various dust models over a range of 
possible compositions and grain size distributions. We also plan to extend our code to handle other plausible 
geometries, such as a dusty disk and a dusty torus. We plan to make these codes publicly available, which 
will provide a useful tool for interpreting available spectropolarimetry results, and in particular, help justify 
new observations that can teach us about the properties of dust in various Galactic and extragalactic 
environments. 

This research was funded by a grant from Israel Science Foundation. We thank the referee John Mathis 
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A. Numerical implementation 

A package of programs was developed to implement the calculations described above. The core of the 
package is the code PRT which calculates the polarized radiative transfer in a plane-parallel dusty slab 
using the adding-doubling method. The light scattering by an optically thick dusty sphere is calculated 
with the code PRT-SPH, using the results, obtained from PRT. The codes were all written in the standard 
ANSI C language, which allows a very efficient use of memory and CPU time. In particular, the memory 
for large multidimensional arrays, like the reflection and transmission matrices, which are a function of the 
number of Gaussian quadrature angles and the number of azimuth modes, is allocated dynamically (i.e. in 
run time). The codes may be ported to other computers where an ANSI C compiler is installed. 



A.l. The code PRT 

The development of PRT was simplified by the availability of the FORTRAN code RT3 generously 



provided by K.F. Evans and G.L. Stephens (Evans & Stephens 1991 ), which solves the polarized radiative 
transfer problem in the plane-parallel case. Some useful ideas and algorithms were adapted from RT3, 
although the realization of our code is completely different from that of RT3 both in the general structure 
and in the details of the code. RT3 was also very valuable for checking PRT during various steps of the its 
development. 

The formulation of the problem that PRT solves is as follows. A plane-parallel slab is illuminated 
by unidirectional radiation from a point-like source and/or diffuse radiation from both above and/or 
below. The slab may be stratified, i.e. composed of a number of layers, each of which may be treated 
as homogeneous, having a specific geometrical thickness and dust grain composition. The dust in each 
layer may contain a number of size-distributed components, which may be spherical grains of various 
types: homogeneous, concentric multilayer and/or composite. The main task of PRT is to compute the 
angular-dependent intensity vectors (Stokes vectors) both for the radiation emergent from the upper and 
lower surface of the slab (external problem) and for the radiation at various levels inside the slab (internal 
problem). PRT solves this problem by the following algorithm: 

1. reads input (see below), and allocates memory for arrays; 

2. calculates quadrature angles and weights; 

3. calculates the scattering matrix for each azimuth mode and each layer; 

4. loop over the azimuth modes: 

(a) loop over the layers: 

i. initializes the doubling by calculating the local reflection and transmission matrices and 
source vectors for an initial very thin sublayer; 

ii. doubling: computes the matrices and sources for the whole layer; 

(b) adding: calculates the matrices and sources for the whole slab; 

(c) calculates the intensity vectors, both external and internal; 

5. converts the intensity vectors from Fourier modes into azimuth space; 
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6. calculates and outputs fluxes of emergent radiation; 

The input includes one main file and a few additional smaller (one per layer) files. A typical main input 
file contains the following parameters: 1. the number of quadrature angles, n^, normally 16-24; 2. the 
number of azimuth modes, n^, normally 4-24, but may be larger if the phase function is strongly forward 
peaked; 3. the number of Stokes parameters, Ust, 1-4; 4. quadrature scheme: Gaussian or double Gaussian; 
5. the flux from the point source, Fq, usually set to 1, and its zenith angle, ^o; 6. the temperature of the 
diffuse radiation, incoming from above and below, usually set to K; 7. the optical thickness of the initial 
doubling layer. Tin, usually 10^^-10"^; 8. the wavelength. A; 9. the number of azimuths, n$, in the range 
0° to 180°, for which output intensity vectors are to be calculated; 10. for each layer: its geometrical depth 
and the name of the additional input file with dust mixture properties. 

A typical additional input file for a layer contains: 1. the dust mass density; 2. number of dust 
components; 3. for each dust component: type of grain constitution (homogeneous, multilayer, composite), 
mass fraction in per cent, name of a file with the size distribution function, the number of grain constituents 
(one for homogeneous grains), and for each grain constituent: name of the constituent and its volume 
fraction in the grain. 

The output includes the angular-dependent intensities and fiuxes of the radiation at various depths 
in the slab. In the simplest case, only the intensities at the top and bottom are printed (refiected and 
transmitted radiation). The results may be stored in text or in binary form. The binary form is especially 
suitable for use by other codes, e.g. when treating the light scattering by an optically thick sphere. 

It is possible to use the Gaussian or the double Gaussian quadrature schemes in PRT. However, we 
have found that in order to produce results which are suitable for usage by PRT-SPH (for calculating light 
scattering by a sphere) the double Gaussian scheme is more preferable, as it generates an angle grid which 
contains very small /i, which is necessary for PRT-SPH to work properly. 

The most important differences between the codes PRT and RT3 are as follow. (1) The memory for 
the multidimensional matrices and vectors in PRT is allocated and freed dynamically, which provides the 
most efficient scheme of memory usage, while in RT3 the arrays are statically allocated. (2) In PRT, the 
scattering matrix for a specified dust grain mixture is calculated using the Mie theory for individual grains 
summing over the mixtures components and grain-size distributions, whereas in RT3 the scattering matrix 
is computed from its expansion in Legendre series in cos 0scat • (3) PRT uses a more efficient realization 
of the FFT procedure with less than 50 lines of C code, whereas the respective FORTRAN code in the 
RT3 contains more than 250 lines. (4) When performing adding and doubling [eqs (^)] and calculating 
the intensities of the radiation inside the slab [eqs (^)], there is a need to work with inverted matrices. 
RT3 directly inverts matrices in these cases. In PRT, we proceed more efficiently. Since the inverted 
matrices in eqs ( ^7| ) and ( ^9|) are combinations like A^^B or A~^b (where A and B are matrices, 6 is a 
vector), it is more computationally suitable to calculate A^^B or A^^b in the following way: first find 
the L[/-decomposition of matrix A and then perform a backsubstitution with the columns of i? or 6 (see 



e.g. Forsythe, Malcolm, & Moler 1977). This is more accurate and saves one matrix multiplication. (5) 
In addition to these differences, PRT was written with the aim of making it more readable and logically 
organized, and thus the general structure of PRT differs from RT3. 

The improved computational scheme of PRT allows standard Mie theory solutions for grains with 
size parameters up to 10^. Thus, we do not need to resort to the Rayleigh-Gans or geometrical optics 
approximations (Bohren & Huffman 1983) for any combination of grain sizes and wavelengths. 
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The accuracy of the resuhs generated by PRT depends on the following parameters: the initial optical 
thickness for doubling Tin, the number of quadrature angles n^, and the number of azimuth modes rim- For 
a proper choice of n^, usually > 16, and n^, the accuracy is set by Tin. For example, setting Tin = 10~^ will 
result in 5-6 digit precision. 

Typical values of running time and virtual memory, required by PRT for various values of ti^ and Um 
are displayed in Table 0. These were derived by using a dusty layer of MRN dust. The results correspond 
to a calculation per one wavelength, A=0.5 /im, and per one angle of incoming radiation, 9q — 30°. The 
optical depth of the layer at A=0.5 fim is ^ 224; Tin for doubling was set to 10~®. The calculations have 
been performed at a DEC Alpha workstation 250''/^^^. The program was compiled with a C compiler GNU 
gcc. 



A.2. The code PRT-SPH 

The code PRT-SPH solves the following problem. A sphere of radius Rs is illuminated by unpolarized 
radiation with flux Fq from a point source (see Figure ||), where the angle between the directions to the 
source and the observer is 8obs- The code PRT-SPH was designed to solve two main tasks: first, the 
computation of the linear and circular polarization maps that requires calculation of the local intensity 
vectors over the surface of the sphere, second, the calculation of the integrated characteristics: the intensity 
and the degree of linear polarization obtained by integrating the local intensities over the sphere. The local 
quantities are calculated using PRT. 

The general algorithm that solves the problem is: 

1. discretize the angles 6 and $ to define small local surface elements; 

2. calculate local intensities and polarization parameters at each of the surface elements (eqs. p3|~|60| and 

I; 



3. calculate the total intensity vector and the degree of linear polarization (eqs. 62-35|). 

The azimuth angle $ is discretized uniformly from —90° to 90° with n$ angles. The zenith angle grid 
is discretized with ng angles in the range 0-180° such that all local elementary squares have equal surface 
area AS'(= 27ri?g/n$ne). The grid angles $fc and 6k are defined by: 



= 90° 



" 2 










, fc = i. 


71$ 





,71$, (Al) 



2fc — 1 

cos6'fc = 1 , k = l,...,no. (A2) 

ne 

The results of the computations for a slab obtained with PRT are output into a binary file, which 
contains the Fourier modes of the intensity vector of the reflected radiation /J'" for a number of fixed grid 
angles 9in and 6'out- This provides the input to PRT-SPH. The local intensity vector ij [in eq. (|5^)] for any 
specific angles ^in, ^out and <i>out is calculated as follows: first, by interpolating the Fourier modes for ^in and 
6'out, next by computing ij with formula (48) for the required <I>out- For the two-dimensional interpolation 



over ^in and Oout, PRT-SPH uses bicubic splines. The accuracy of this approach for calculating l| generally 
lies well within 0.1% when a grid 16x16 for fiin and /lout is used for the PRT data. 
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To compute the integral parameters Igph and Qsph with formulae (|6^|6^), the respective integrals are 
approximated by sums over the surface elements. For the grid of 40x90 for 6 and $ over the sphere, typical 
in our calculations, the accuracy of the results is better than 0.01-0.1%. 



B. Code Testing 

B.l. Comparison with analytic results for Rayleigh scattering 

To test the accuracy of PRT and PRT-SPH, we performed calculations for a pure Rayleigh scattering 
sphere. This corresponds physically to a Thomson scattering electrons sphere. This case was carefully 



studied and there are various analytical results to compare with, e.g. Chandrasekhar (1960) 



We calculated the inclination dependence of the intensity and the degree of linear polarization for 
radiation transmitted through an optically thick slab. We used n^=36, n,„=4, rin=10~^ and t=1460 for 
the total optical depth of the scattering layer. The calculations made with PRT agree with the analytical 



results of Chandrasekhar (1960) to better than 0.05% 



The optical depth dependence of the flux of radiation transmitted through a Rayleigh scattering slab 
is shown in Figure |l8|, as calculated with PRT. For comparison, we plot the large-r asymptotic analytic 
solution given by Yanovitskij (1997): 

We find an agreement to within 1% between the analytic and numerical solutions at r >3. 

As a side note, we tried to use our numerical solution, which is exact at all r, to obtain an improved 



analytic approximation which generalizes equation Bl to a larger range of r. We found a surprisingly good 
approximation to our result which is valid for all t. 

, , 4 1.265 -a 

Fa T = (B2) 

' 3r+ 1.423-6 ^ ' 

with a=1.019 and 6=1.194. The accuracy of this simple approximation is 0.3-2% for r < 0.5, 2-4% for 
0.5 < r < 5, and 0.5-2% for 5 < r. 

In a pure scattering atmosphere no flux is lost due to absorption, and thus energy conservation implies 
that the incident flux must equal the sum of transmitted and scattered fluxes. The right panel in Figure 
shows the accuracy to which energy is conserved in PRT. There error increases roughly linearly with r 
due to accumulation of errors in the doubling procedure. However, the errors for t ^ 100 — 1000, typical 
for this study, are stiU very smaU (<0. 0001-0.001%). 



B.2. Comparison with previous results for a dusty sphere 



White ( 1979b| ) calculated the scattering properties of an optically thick dusty sphere with an MRN 
composition, also using a doubling method for the radiative transfer. Figure ^ shows a comparison of the 
wavelength-dependent spherical albedo, maximum linear polarization, and ratio of forward to backward 
scattering, taken from the figures in White (1979b) to those calculated with PRT and PRT-SPH. The 
overall consistency is quite satisfactory, considering the significant improvements in the optical constants 
for graphite and silicate grains since the study of White (1979b). 
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C. Electron scattering versus dust scattering sphere 

Figure ^ presents maps of the linear polarization of an optically thick electron scattering sphere, 
obtained with PRT-SPH. Comparison with corresponding maps of dusty spheres (Figs |l3|-|l^ indicates 
that, as expected, multiple scattering has a larger contribution to the polarization in the electron scattering 
sphere. For example, at 8obs=180°, the local degree of p reaches 7-8%, which is twice the value for a 
dusty sphere, and a similar effect is seen at very small Qohs- On the other hand, the maximum local 
polarization (at 8obs=90°), is 55-60% which is significantly lower than the maximum values of up to 90% 
(for A=0.05 /im) found for a dusty sphere. Dust absorption enhances the polarization since it suppresses 
the contribution of multiple scattering, which adds unpolarized light, or light which is polarized at right 
angle to the single scattered light (for 8obs=90°), in both cases lowering the total polarization. 



Figure 21 presents the integrated scattering characteristics of the electron scattering sphere in 
comparison to those of a dusty sphere. The electron scattering sphere scatters light more efficiently at most 
angles (8obs >20°). Its albedo is essentially unity, versus 4 — 11% only for a dusty sphere. However, a dusty 
sphere is a much more efficient polarizer, allowing maximum Psph of 70-80%, versus only about 30% for the 
electron scattering sphere, again due to the greater contribution of multiple scattering in the later case. 

Note that the integrated electron scattering sphere polarization is rotated by 90° for 8obs=0-20°. This 
effect is not seen in the Monte Carlo calculations of Code & Whitney (1995), most likely due to the very low 
intensity of light scattered at this range, coupled with the low polarization amplitude. This demonstrates 
the advantage of an accurate solution over Monte Carlo calculations, though the former is limited to the 
T » 1 case only. 
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Table 1. CPU time and virtual memory, required by PRT. 
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Fig. 1 . — Schematic representation of light scattering by a slab of optical thickness tq . Unpolarized radiation 
with a flux _Fo is incident from a zenith angle 6q. The downward and upward intensities (/^ and J^) are a 
function of t, /j, — cos9 and azimuth angle Diffuse radiation may be incident from the top (/^ ) and the 
bottom (7j[). The radiation which emerges from the top (ij) and bottom (/^), results from the transfer which 
is expressed by the reflection (i?^^, R^^) and transmission {T^^, i?TT) matrices, plus a possible contribution 
of internal sources {S^, S^). 

Fig. 2. — The scattering geometry for a plane-parallel slab. The slab is in the XY plane and Z is the 
normal. The incoming beam makes an angle with the Z axis, and the outgoing beam is in the direction 
^out , ^out ; where #out is the azimuth difference between the incoming and outgoing beams. The scattering 
angle Oscut is the angle between the continuation of the incoming beam A'"0 and the outgoing beam OA°^^ . 
The rotations angles of the incoming (outgoing) beam ii (22) is measured from the meridian plane of the 
incoming beam A^^OZ (the plane A°^'^OZ), specified by nJJJ (n™'), and the scattering plane A'"OA°"*, 
specified by n™ (n°"*). Right panel: A view of the scattering geometry along the outgoing beam, showing 
the plane of linear polarization of single scattered radiation, and of multiply scattered radiation. The first 
is generally at right angle to the plane defined by the incoming and outgoing beams, while the second is 
generally at right angle to the plane of the slab. 

Fig. 3. — The scattering geometry for a dusty sphere. The direction of the incoming beam is A^'^S, and 
it gets scattered by Gobs along the direction 5^4°"*. The scattering point S coordinates on the sphere are 
zenith angle 6 and azimuthal angle $ (measured from the X axis), and the local normal is along n{6,^). 
The scattering at S is calculated using the local slab approximation, as presented in Figure 2. 

Fig. 4. — The single scattering albedo and the phase function asymmetry parameter g as a, function of 
wavelength. Note the sharp drop in albedo below 2200 A, and the tendency of scattering to become strongly 
forward peaked at short wavelengths. 

Fig. 5. — The extinction curve for the MRN dust mixture compared to the mean Galactic extinction curve 
for Rv=3.1 (Cardelli, Clayton, & Mathis 1989). Note that the Galactic extinction at 1000 A is somewhat 
uncertain. 

Fig. 6. — The scattering and polarization properties of optically thin MRN dust. Left panels show the 
A dependence of the normalized intensity /(A, Gscat)/-'^(1 /um,90°), degree of linear polarization, p{X), and 

polarized intensity I x p. Right panels show the ©scat dependence for the same parameters. The polarization 
plane is generally at right angle to the scattering plane, excluding the cases where p(A) < (marked by bold 
lines and filled points in the lower panels), which refer to polarization parallel to the scattering plane. 

Fig. 7. — The logarithm of the total intensity /t(A) of light scattered by an optically thick MRN dusty slab. 
Each row corresponds to a given Hin, each column to a given ^out, and each subpanel presents curves for a 
range of Mout- The incident flux is unity. 

Fig. 8. — The degree of linear polarization for the same parameters as in Fig. 7. 

Fig. 9. — The position angle of linear polarization for the same parameters as in Fig. 7. 

Fig. 10. — The degree of circular polarization for the same parameters as in Fig. 7. Note that due to 
symmetry the circular polarization is zero at $out=0° and at 180°. 

Fig. 11. — The linear polarization of the light backscattered (i.e. /Xin = Mout) ^out = 180°) by an optically 
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thick slab as a function of A and /i. Only multiple scatterings contribute to p in this case, and its amplitude 
is generally smaller by a factor of 10-20 compared with the single scattering contribution. 

Fig. 12. The plane albedo of an optically thick dusty MRN slab as a function of A and fi-m- The albedo 
of an optically thick sphere is also indicated. 

Fig. 13. — Maps of the linear and circular polarization of light scattered by an optically thick MRN dusty 
sphere at A = 0.05 /zm for various Qobs- Note the very large differences in the relation between arrow size 
and polarization amplitude at different 6obs- 

Fig. 14.- As in Figure 13 for A = 0.2124 /xm. 

Fig. 15. — As in Figure 13 for A = 1.0 fim. 

Fig. 16. The scattering and polarization propertic;s of an optically thick MRN dusty sphere. All parameters 
are as shown in Fig. 6, excluding the intensity which is not normalized here. 

Fig. 17. — Comparison of scattering by optically thin dust versus an optically thick dusty sphere. Upper 
two panels show the A dependence, and the lower two panels the 6 dependence. Note the large difference in 
the wavelength dependence of 7(A, 0)/J(l /xm, 6) for the optically thin and thick cases, but the remarkable 
similarity of p(A). 

Fig. 18. — The flux of the radiation transmitted through a Rayleigh scattering slab as a function of optical 
thickness {left panel) as calculated with PRT. The incident flux is unity with /iin=l. The large-r analytic 
asymptotic solution j^=| ^^'^423 of Yanovitskij (1997), agrees well with our results. We also show our 
improved approximation F=^:f^^§^ with a=1.019 and 6=1.194, which is accurate to 4.0% for all r. The 
accumulated error in flux conservation, as calculated with PRT for the Rayleigh scattering atmosphere, is 
shown in the right panel. 

Fig. 19. — Comparison with the earlier calculations of White (1979b) for the spherical albedo, maximum 

linear polarization, Pmaxj and ratio of forward to backward scattering, p, for scattering by an optically thick 
dusty sphere of the MRN dust composition. The agreement is acceptable, given the updates in the dielectric 
function since the study of White. 

Fig. 20. — Maps of linear polarization of the light scattered by an optically thick Rayleigh scattering sphere 
(electron scattering) for various ©obs- 

Fig. 21. — Comparison of the scattering properties of an optically thick electron scattering sphere versus an 
optically thick dusty sphere. The intensity of scattered light is generally much higher for an electron sphere, 
but its integrated polarization is significantly lower. 
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